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Abstract 

We develop sub-Nyquist sampling systems for analog signals comprised of several, possibly overlapping, finite 
duration pulses with unknown shapes and time positions. Efficient sampling schemes when either the pulse shape 
or the locations of the pulses are known have been previously developed. To the best of our knowledge, stable and 
low-rate sampling strategies for continuous signals that are superpositions of unknown pulses without knowledge 
of the pulse locations have not been derived. The goal in this paper is to fill this gap. We propose a multichannel 
scheme based on Gabor frames that exploits the sparsity of signals in time and enables sampling multipulse signals 
■ at sub-Nyquist rates. Moreover, if the signal is additionally essentially multiband, then the sampling scheme can be 

adapted to lower the sampling rate without knowing the band locations. We show that, with proper preprocessing, 
the necessary Gabor coefficients, can be recovered from the samples using standard methods of compressed sensing. 
In addition, we provide error estimates on the reconstruction and analyze the proposed architecture in the presence 
^ ■ 

i of noise. 
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I. Introduction 
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^C) , One of the common assumptions in sampling theory suggests that in order to perfectly reconstruct a bandlimited 



analog signal from its samples, it must be sampled at the Nyquist rate, that is twice its highest frequency. In practice, 
however, all real life signals are necessarily of finite duration, and consequently cannot be perfectly bandlimited, 
due to the uncertainty principle Q). The Nyquist rate is therefore dictated by the essential bandwidth of the signal, 
that is by the desired accuracy of the approximation: the higher the rate, meaning the more samples are taken, the 
better the reconstruction. 

In this paper we are interested in sampling a special class of time limited signals: signals consisting of a stream 



of short pulses, referred to as multipulse signals. Since the pulses occupy only a small portion of the signal support, 
intuitively less samples, then those dictated by the essential bandwidth, should suffice to reconstruct the signal. 

There are two standard approaches in the literature to sample such functions. One is to acquire pointwise samples 
and approximate the signal using Shannon's interpolation formula J5), 0. The reconstruction error can be made 
sufficiently small with just a finite number of samples, when the signal is sampled dense enough. However, this 
strategy results in many pointwise samples that are zero, leading to unnecessary high rates. The second, is to collect 
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Fourier samples and approximate the signal using a truncated Fourier series. However, the Fourier transform does 
not account for local properties of the signal, hence this method cannot be used to exploit signal structure and 
reduce the sampling rate. Both strategies require the Fourier transform of the signal to be integrable and do not take 
the sparsity of the signal in time into account. Moreover, exact pointwise samples needed for Shannon's method 
requires implementing a very high bandwidth sampling filter. Here we show that these problems can be alleviated 
using Gabor frames H. 

Gabor samples, which are inner products of a function with shifted and modulated versions of a chosen window, 
are a good compromise between exact pointwise samples and Fourier samples. In particular, we show that all square- 
integrable time limited signals, without additional conditions on their Fourier transforms, can be well approximated 
by truncated Gabor series. Furthermore, Gabor samples, taken with respect to a window that is well localized in time 
and frequency, provide information about local behavior of any square integrable function and reflect the sparsity 
of a function either in time or frequency. The price to pay is a slightly greater number of samples necessary for 
approximation, that comes with using frames, namely, overcomplete dictionaries. The use of frames is a result of 
the fact that Gabor bases are not well localized in both time and frequency [5 |. In all three approaches (pointwise, 
Fourier, Gabor) the number of samples necessary to represent an arbitrary time limited signal is dictated by the 
essential bandwidth of the signal and the desired approximation accuracy. 

Recently, there has been growing interest in efficient sampling of multipulse signals fS), Q, O, (9J. This 
interest is motivated by a variety of different applications such as digital processing of certain radar signals, which 
are superpositions of shifted and modulated versions of a single pulse Q, ifToll . ifTTIl . Another example is ultrasound 
signals, that can be modeled by superpositions of shifted versions of a given pulse shape [9]. Multipulse signals 
are also prevalent in communication channels, bio-imaging, and digital processing of neuronal signals. Since the 
pulses occupy only a small portion of the signal support, intuitively less samples should suffice to reconstruct the 
signal. 

Prior works mentioned above assumed that the signal is composed of shifts of a single known pulse. Such signals 
are completely characterized by a finite number of parameters and fall under the class of finite rate of innovation 
(FRI) signals introduced in [6 |. The sampling schemes proposed in [8] operate at the minimal sampling rate required 
for such signals, determined by the rate of innovation |6|. In this case without noise, perfect recovery is possible 
due to the finite dimensionality of the problem. 

In this paper we consider sampling of multipulse signals when neither the pulses nor their locations are known. 
The pulses can have arbitrary shapes and positions, and may overlap. The only knowledge we assume is that our 
signal is comprised of N pulses, each of maximal width W. Despite the complete lack of knowledge on the signal 
shape, we show that using Gabor frames and appropriate processing, such signals can be sampled in an efficient and 
robust way, using far fewer samples than that dictated by the Nyquist rate. The number of samples is proportional to 
WN, that is, the actual time occupancy. More precisely, we need about 4// _1 f2'W A A< r samples, where fl' is related 
to the essential bandwidth of the signal and // G (0, 1) is the redundancy of the Gabor frame used for processing. 
When the signal is additionally sparse in frequency with only S essential bands of width no more than flw, the 
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sampling rate can be further reduced. For such signals, we need about Spi^^^ft'yyWNS samples, where U! w is 
related to the width Qyj/ of the essential bands of the signal. In contrast, Nyquist-rate sampling in both settings 
requires about 0'/3 samples, where (3 is the signal duration. If the signal occupies only a small portion of its time 
duration, such that 4/j,~ 1 WN <C /?, respectively 2il' w S <C f2, then our scheme results in a substantial gain over 
Nyquist-rate sampling. 

The sampling criteria we consider are: a) minimal sampling rate that allows almost perfect reconstruction, b) no 
prior knowledge on the locations or shapes of the pulses, and c) numerical stability in the presence of mismodeling 
and noise. To achieve these goals we combine the well established theory of Gabor frames |4) with compressed 
sensing (CS) methods for multiple measurement systems fl2l . fl3l . 03). Our scheme consists of a multichannel 
system that modulates the input signal in each channel with a parametric waveform, based on a chosen Gabor 
frame, and integrates the result over a finite time interval. We show that by a proper selection of the waveform 
parameters, the Gabor samples can be recovered, from which the signal is reconstructed. We also consider the case 
in which the signal exhibits additional sparsity in frequency, as is common in radar signals, and show that using 
our general scheme the sampling rate can be further reduced. To recover the signal in this case we solve two CS 
problems. We then prove that the proposed system is robust to noise and model errors, in contrast with techniques 
based on exact pointwise samples. 

Our development follows the philosophy of recent work in analog CS, termed Xampling, which provides a 
framework for incorporating and exploiting structure in analog signals to reduce sampling rates, without the need for 
discretization lfl5l . fl6l . Xampling combines standard analog sampling methods with CS digital recovery techniques. 
A pioneer sub-Nyquist system of this type is the modulated wideband converter (MWC) introduced in ifTTIl based 
on the earlier work of ifTSl . This scheme targets low rate sampling of multiband signals. Sub-Nyquist sampling is 
achieved by applying modulation waveforms to the analog input prior to uniformly sampling at the low rate. 

Another system that falls into the Xampling paradigm is that of [ 8 1 which treats multipulse signals with a known 
pulse shape. The proposed sampling scheme is based on modulation waveforms as in the MWC. However, while in 
the MWC the modulations are used to reduce the sampling rate relative to the Nyquist rate, in [ 8 1 the modulations 
serve to simplify the hardware and improve robustness. 

Gabor frames were recently used to sample short discrete pulses in fl9l . The authors analyzed standard CS 
techniques for redundant dictionaries, and applied their results to radar-like signals. Finite discrete multipulse signals 
were also treated in [20] where the authors modeled the signals as convolutions of a sparse signal with a sparse 
filter, both sparse in the standard basis of C N . The important difference between |[T9l , ll20l and our work is that 
the former handles discrete time signals. In contrast, our method directly reduces the sampling rate of continuous 
time input signals without the need for discretization. 

The paper is organized as follows. In Section HI] we introduce the notation and basic problem definition. Since 
the main tool in our analysis is Gabor frames, in Section [HI] we recall basic facts and definitions from Gabor 
theory and show that truncated Gabor series provide a good approximation for time limited functions. Based on 
this observation, in Section Hy] we introduce a sub-Nyquist sampling scheme for multipulse signals. In Section M 
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we show that our system can also be used to efficiently sample radar-like signals, who are sparse in time and 
frequency. Section [VT] points out connections to recently developed sampling methods, while Section [VHl is devoted 
to implementation issues. The important part of our design are Gabor windows, which we review in Section [Villi 
In particular, we summarize several methods to generate compactly supported Gabor frames. We demonstrate our 
theory by several numerical examples in Section [IX] 

II. Problem Formulation and Main Results 

A. Notation 

We will be working throughout the paper with the Hilbert space of complex square integrable functions L,2(M,), 
with inner product 

f(t)gjtjdt for all /,jeL 2 (R) 

-oo 

where g(t) denotes the complex conjugate of g(t). The norm induced by this inner product is given by ||/||| = (/ 5 /). 
The Fourier transform of / 6 is defined as 

/oo 
f(t)e- 2niut dt 
-00 

and is also square integrable with \\ f\\-2 = l/lh- 

A main tool in our derivations are Gabor frames, which we review in Section IIII-AI Two important operators 
that play a central role in Gabor theory, are the translation and modulation operators defined for i,weias 

T x f(t) := f(t - x) , M u f(t) := e 2 ^ l f(t) , 

respectively. The composition M u T x f(t) — e 2vlult f(t ~ x) is called a time-frequency shift operator and gives rise 
to the short-time Fourier transform. For a fixed window g £ L2CR), the short time Fourier transform of / e L2GR) 
with respect to g(t) is defined as 

V g f(x,u) := (f,M u T x g). 

Many derivations, and especially input-output relations for our sampling systems, will be presented in the compact 
form of matrix multiplications. We denote matrices by boldface capital letters, for example C, D, and vectors by 
boldface lower case letters, such as x, z. 

Our recovery method relies on CS algorithms. An important notion in this context is that of the restricted isometry 
property (RIP). A matrix C is said to have the RIP of order S, if there exists < 5 < 1 such that 

(l-5)||u||I<||Cu||i<(l + J)||u||| 

for all 5— sparse vectors u ET1 . 
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Fig. 1: Schematic example of a multipulse signal with N — 6 pulses each of width no more than W. In the example, 
two of the pulses are overlapping. 



B. Problem Formulation 

We consider the problem of sampling and reconstructing signals comprised of a sum of short, finite duration 
pulses. A schematic representation of such a signal is depicted in Fig. Q] We do not assume any knowledge of the 
signal besides the maximum width (support) of the pulses. More formally, we consider real valued signals fit) of 
the form 

N 

f(t)=y h n (t) , where max Isupp h n \ < W . (1) 

n=l 

The number of pulses N and their maximal width W are assumed known. The pulses may overlap in time, as in 
Fig. Q] We assume that /(f) is supported on an interval [— /3/2,/3/2] with NW <C (3. Our goal is to recover fit") 
from the minimal number of samples possible. 

Due to the uncertainty principle, finite duration functions cannot be perfectly bandlimited. However, in practice the 
main frequency content is typically confined to a finite interval. We refer to such signals as essentially bandlimited. 
More formally, we say that fit") is essentially bandlimited, or en— bandlimited to F = [—0/2,17/2], if for some 
en < 1 

\ 1/2 

\J{u:)\ 2 du) <en||/|| 2 . (2) 



The symbol F c denotes the complement of the set F. The adjective 'essential' refers to the fact that the energy of 
f{uS) outside [ — 0/2, Sl/2] is very small. We denote the set of multipulse signals ([TJ timelimited to [— B/2,/3/2] 
and essentially bandlimited to [-0/2, 0/2] by MV(N, W, j3, O). 

There are three interesting special cases that fall into the model (Q~|i. The first is when h n (t) are shifts of a known 
pulse hit), so that h n (t) = a n h{t — t n ) for some t n , a n € R. In this case, the problem is to find 2N parameters, 
the amplitudes <t„ and shifts t n . This setting can be treated within the class of finite rate of innovation problems 
[6 1, [8 1, |9|. We return to this scenario in Section [VT1 and discuss the relation to our work in more detail. A second 
class, is when the location of the pulses h n (t) are known but the pulses themselves are not. The third, most difficult 
scenario, is when neither the locations nor the pulses are known. Our goal is to develop an efficient, robust, and 
low-rate sampling scheme for this most general scenario. We will later see that our system can be used to sample 
signals from the other two cases as well, at their respective minimal rates. In Section [V] we show that our system 
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can be additionally used to reduce the sampling rate of a special subclass of AiViN, W, (3, 51), which are multipulse 
signals whose frequency content is concentrated on only a few bands within [ — 0/2, 0/2]. 

We aim at designing a sampling system for signals from the model JAV(N, W, f3, £1) that satisfies the following 
properties: 

(i) the system has no prior knowledge on the locations or shapes of the pulses; 

(ii) the number of samples should be as low as possible; 

(iii) the reconstruction from the samples should be simple; 

(iv) the original and reconstructed signals should be close. 

C. Main Results 

The proposed multichannel sampling method, depicted in Fig. @] is a mixture of ideas from Gabor theory and 
Xampling [16|. It consists of a set of modulators with functions qj yTn (t), followed by integrators over the interval 
[— /3/2, /3/2]. The system depends on an appropriately chosen Gabor frame with redundancy degree [i G (0,1), 
generated by a compactly supported window that is well localized in the frequency domain. This frame provides a 
sparse representation for Ai"P(N, W, (3, Q). The modulating waveforms qj. m (t), formally defined in ©, are different, 
finite superpositions of shifted versions of the chosen Gabor window. The goal of the modulators is to mix together 
all windowed pieces of the signal with different weights, so that, a sufficiently large number of mixtures will allow 
to almost perfectly recover relatively sparse multipulse signals. The resulting samples are weighted superpositions 
of Gabor coefficients of the signal with respect to the chosen frame. CS methods fl2l . lfl3l . lfT4l are then used to 
recover the relevant nonzero signal coefficients from the given samples. 

The number of rows in the resulting CS system is about 4iV/x _1 ; it is a function of the number of pulses 
present in the signal and the redundancy /i of the frame. Since CS algorithms are used to recover the relevant 
coefficients, the exact number of rows is dictated by the RIP constant of the matrix containing the coefficients of 
the waveforms, and is given by 0{AN yT 1 log(/3' / (4NW))). In the case of purely multipulse signals, the number 
of columns is a function of the desired accuracy of the approximation, and equals about Vl'W . However, when 
the signal is essentially mutiband, with S bands of width VLw, then the number of columns can be reduced to 
about 2il' w WS, proportional to the actual frequency content of the signal. Again, since CS methods are used in 
the recovery process, the overall number of columns is dictated by the RIP constant of the matrix containing the 
coefficients of the waveforms, and is given by 0(2fl' w WS \og(fl' w / (2Q' W S))) . The quantities /?', SI' and Q! w are 
related to (3, SI and £l' w , respectively, and depend on the chosen Gabor frame. 

After finding the Gabor coefficients, we recover the signal using a dual Gabor frame. The function f(t) recon- 
structed from the post-processed coefficients satisfies 

||/- fh < &o(en + e B )\\fh + ^iIMb + C 2 ||n 2 || a , 

where f(t) is the original signal, Co is a constant depending on the Gabor frame, and es is related to the essential 
bandwidth of the chosen Gabor window. The first term is due to the signal energy outside the essential bandwidth. 
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The values of ni and n 2 reflect the noise level in the signal (mismodeling error) and the samples, respectively, 
while the constants C\ and C 2 depend on the CS method used for recovery of the Gabor coefficients. If f(t) is 
perfectly multipulse and the sampling system is noise free, then ni = n 2 = 0. For multipulse signals that are 
essentially multiband, ni is related to the signal energy outside the essential bands of f(t). 

III. Sampling using Gabor frames 

We begin by recalling some basic facts and notions from Gabor theory that will be used throughout the paper, and 
then show how Gabor frames can be used to sample multipulse signals with known pulse locations. In Section HVl 
we expand the ideas to treat the unknown setting. 

A. Basic Gabor Theory 

A coUection Q(g, a, b) = {M u T ak g(t) = e 27rMt g(t -ah); k,l £ 1} is a Gabor frame for L 2 (M) if there exist 
constants < A\ < A 2 < 00 such that 

^i||/|| 2 < J2^f' M ^ T ^9)\ 2 <A 2 \\f\\ 2 

for all / £ L 2 (R). The frame is called tight, if Ai = A 2 . By simple normalization every tight frame can be changed 
to a tight frame with frame bounds equal to one. Therefore, when we talk about tight frames we will mean frames 
with frame bounds A\ = A 2 = 1. Every signal / e L 2 (M.) can be represented in some Gabor frame [4]. 

A Gabor representation of a signal /(t) comprises the set of coefficients {zk.i}k,iei obtained by inner products 
with the elements of some Gabor system Q(g, a, b) 0J: 

*k,l = (/, M bl T ak g) = e 27rafcw (/, M_ ak T bl g) . 

The coefficients z^; are simply samples of a short-time Fourier transform of f(t) with respect to g(t) at points 
(ak, bl). If G{g, a, b) constitutes a frame for L 2 (R), then there exists a function 7 £ L 2 (M) such that any / £ L 2 (M) 
can be reconstructed from {^,i}fe,iez using the formula 

f=J2 Zk,lM bl T akl . (3) 
fc.iez 

The Gabor system ^(7,0,6) is the dual frame to Q(g,a,b). Consequently, the window j(t) is referred to as the 
dual of g(t). Generally, there is more than one dual window -f(t). The canonical dual is given by 7 = S~ x g, where 
S is the frame operator associated with g(t), and is defined by 5/ = ^ klt - z (f 1 M b iT a kg)MuT a kg. There are 
several ways of finding an inverse of S, including the Janssen representation of £, the Zak transform method or 
iteratively using one of several available efficient algorithms Bl . 

Here we will only be working with Gabor frames whose windows are compactly supported on some interval 
[—a/2, a/2] and lattice parameters a = fia, b = 1/a for some fi £ (0, 1). For such frames, the frame operator 
takes on the particularly simple form 
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Fig. 2: The relation between / and the shifts of the support of g when /i = 1. When sappTwkd for some k, does 
not overlap any of the pulses of /, then z^.i = for all I. 

The frame constants can be computed as A\ = ess inf S(t) and = ess sup S(t). The canonical dual is then 
7(i) = 6S' _1 (i)g(t). For tight frames the dual atom is simply 7(i) = A^ 1 bg(t). A necessary condition for Q(g, a, b) 
to be a frame for ^(R) is that a b < 1, while Gabor Riesz bases can only exists if ab = 1 |4|. Thus the ratio 
l/(ab) measures the redundancy of Gabor systems. 

Since one key motivation for considering Gabor frames is to obtain a joint time-frequency representation of 
functions one usually attempts to choose the window g(t) to be well localized in time and frequency. While the 
Balian-Low theorem [5] makes it impossible to design Gabor Riesz bases with good time-frequency localization, 
it is not difficult to design Gabor frames with excellent localization properties. For instance, if g(t) is a Gaussian, 
then we obtain a Gabor frame whenever ab < 1. Therefore, to obtain a well localized window one needs to allow 
for certain redundancy. In Section [VIIII we discuss in detail how to construct frames and their duals with compactly 
supported windows based on |22|, |23l . 

We consider windows g(t) that are members of so-called Feichtinger algebra, denoted by So \24\. Such windows 
guarantee that the synthesis and analysis mappings are bounded and consequently result in stable reconstructions, 
and that the dual window is in Sq. Formally, 



where ip(t) — e _7ri . The norm in So is defined as ||/||s := ||VL,/||i. Examples of functions in So are the Gaussian, 
B-splines of positive order, raised cosine, and any L\ (R) function that is bandlimited or any (R) function that is 
compactly supported in time with Fourier transform in Xi(R). Note, that the rectangular window is not a member 
of So since its Fourier transform is not in ii(K). 

B. Truncated Gabor Series 

It is well known that time limited L2(R) functions, whose Fourier transform is additionally in Li(R), can be 
well approximated with a finite number of samples using a Fourier series. We now show that the same is true for 
Gabor series, without assuming anything additional on the signal besides that it is square integrable. 

Let Q(g, a, b) be a Gabor frame with g(t) compactly supported on an interval [—a/2, a/2], a = [ia and b = 1/a 
for some /j£ (0,1). The reason for using compactly supported windows is that for every function f(t) time limited 
to [— (3/2, /3/2], the decomposition of (f3]l reduces to 





(4) 
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where j(t) is a dual window and Kq denotes the smallest integer such that the sum in contains all possible 
non-zero coefficients Zk : i- The exact value of Kq is calculated by 

~P + a' 



~ + (K + l)a> |=>«b = 



2a 



1 . 



The number of frequency samples I necessary for e reconstruction is dictated by the pair of dual windows (<?, 7) 
as incorporated in the following theorem, which is an extension of Theorem 3.6.15 in l24l . 

Theorem III.l. Let f(t) be a finite duration signal supported on the interval [— /3/2,/3/2] and en— bandlimited 
to [—fl/2,Q/2] and let Q(g,a,b) be a Gabor frame described above with the dual atom 7 £ So- Then for every 
tB > there exists an Lq < 00, depending on the dual window j(t) and the essential bandwidths of g(t) and 
f(t), such that 



K L 

f~ 2J X! z kd M blT a kl 
k— — Kq I— — L n 



<C (e n + e s )||/|| 2: 



where Cq = blMlsollfl'II.S'o w ^ tn C a ,b = (1 + l/a) 1 / 2 (l + 1/b) 1 / 2 a constant depending on the chosen Gabor 
frame. 



Similar estimates also appear in [25]. 

Proof: See Appendix lAl ■ 
The exact number of frequency coefficients Lq is dictated by the essential bandwidth of g(t). More precisely, 
if g c (t) is a [— B/2, B/2}— bandlimited approximation of g(t) in So, that is \\g — g c \\s ^ e s||9lls ' tn en Lq = 

Theorem IIII. 1 1 states that finite duration, essentially bandlimited signals, can be well approximated using just the 
dominant coefficients in the Gabor representation. 

The number of samples depend on the chosen frame and the accuracy of the approximation. To minimize the 
number of samples for a chosen accuracy of the approximation, we select p > 1/2 (which reduces the number of 
samples in time) and construct a window that is well localized in frequency (which reduces the number of samples 
in frequency). Therefore, there is an interplay between the number of samples in the frequency domain and the 
number of samples with respect to time. The total number of Gabor coefficients, meaning samples of the short-time 
Fourier transform, is related to a somewhat larger interval [— /3'/2, /3'/2] C [— j3/2, (3/2], where j3' — (3 + a, with 



K = 2K + 1 
with L = 2Lq - 



—, in the time domain and a larger interval [-&,' /2, Q'/2] C [-fi/2, tt/2], where Q' = + B, 
1 w t~, in the frequency domain. Overall, the required number of samples is 



KL 





"/3 + a" 




~n + B~ 










-0 




2a 


2b 



a b 

When p is close to one, and g(t) is well localized in frequency forming a tight frame, the number of required 
samples is close to 17/3. For a fixed /1 and a chosen accuracy of approximation, the number of frequency samples in 
a tight frame depends on the decay properties of Therefore, to minimize the number of channels, we need to 
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Fig. 3: Schematic representation of a matrix Z of Gabor coefficients in the case of multipulse signals, on the left, 
and multipulse essentially multiband signals, on the right. Empty circles denote zero values and crosses denote 
small, but nonzero values of Zk,i> 
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choose a window g(t) that exhibits good frequency localization. On the other hand, having already chosen a frame 
Q(g,a,b), if we desire to improve the accuracy of approximation, then the number Lq of 'frequency' coefficients 
has to increase. 

C. Multipulse Signals with Known Pulse Locations 

If a <C /?, and the signal is multipulse, then many of the K Gabor coefficients are zero. Indeed, if the shift 
git — ak) does not overlap any pulse of f(t) then 

rP/z 

Zk,l = (/, M bl T ak g) = / f(t)g(t - ak) e - 2 * Mt dt = , 

J -fill 

for all / € Z. Therefore, when the locations of the pulses are known, we can reduce the number of samples from 
KL to ML, where M < K is the number of ks, \k\ < Kq, for which Zk.i 7^ 0. To reduce M to minimum, one 
needs to choose a Gabor frame that allows for the sparsest representation of /(£) with respect to the index k. 

For signals from AiV{N, W, (3, Q), an optimal choice is an atom g(t) that is supported on [— W/2,W/2] and 
shift parameters a = fiW, b = 1/W for some fj, E (0, 1). In that case at most [2/i -1 ] shifts of g(t) by ak — pWk 
overlap one pulse of f(t). Indeed, when \i = 1, at most two shifts of g(t) overlap one pulse, as depicted in Fig. [2] 
When fi < 1, at most shifts of supp g overlap one pulse of f(t). This can be calculated from 

^ < + nWK x Ki > 

2 2 K x - K 2 > 2/i- 1 . 

-f > f + pWK 2 ^K 2 < -n- 1 

Let Z denote the K x L matrix of dominant Gabor coefficients. For functions / 6 .M"P(iV, W, (3, f2) each column 

Z[l] = [z-k ,i, ■ ■ ■ , zk ,i] T of Z has at most [2/x" 1 ] nonzero entries. Moreover, all columns Z[l] have nonzero 

entries at the same places, as modulations e 27rlblt applied to f(t) do not change the positions of the pulses. The 

matrix Z is schematically depicted in Fig. [3] Therefore, the necessary Gabor coefficients can be obtained with only 

ML « 2VL'WNii- 1 channels, where M = \2^]N and L w Vt'W. 

D. Method Comparison 

Since time limited functions can be reconstructed only to a certain accuracy, we refer to the minimal number 
of samples as the minimal number required to reconstruct the signal with a desired accuracy. For an e accuracy 
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of approximation using the Fourier series and Shannon's interpolation Q, @ methods, the minimal number of 
samples is of order Q"/J, where f2" is such that 

j/(w)|dw< ell/Ha, 

with Fx = [— f2"/2, Q"/2] and we have to assume that / £ Li(R). For the above to be satisfied with e = en, as in 
(O, ft" has to be greater than il since H/H2 < II /111- The approximation error using Fourier series is then given by 

/(*)- E f (T) e2rft/ 1 2 ^ E ift/fli < 5 > 

|Z|<£ ^ ' |J|>i 

where / are the Fourier coefficients and L has to be equal at least il"/3/2 to achieve e approximation. When 
the signal is multipulse, L cannot be reduced because the Fourier transform does not account for local signal 
properties. 

The approximation error using Shannon's interpolation formula equals 

\f(t) -(Sn» /)(*)! < / l/MI du, (6) 

J\u\>Q"/2 

where Kq is the largest integer less then Q,"j3/2 and 

(&>«/)(*) = E /(^)sinc(^"(i-fc/0")). 

For fe > Kq, f (-St) — as f(t) is of finite duration, so that about Q,"/3 pointwise values of f(t) must be evaluated 
to achieve e accuracy. If the signal is multipulse and the pulse locations are known, then this number can be reduced 
to NWVt" samples, with WQ" samples per pulse. 

For a Gabor frame with redundancy /1, we achieve e approximation with a minimal number of samples of order 
51'/3'/i _1 as long as the Gabor window g(t) and its dual j(t) are such that 

/ \ l/ 2 

/ / \V g f{x,w)\ 2 dxdu) < I \\fh, 
\JeJfz J ^aMnWSo 

where E = [-ft /2, /3' /2] and F 2 = /2,ti' /2}. The Q' is an enlargement of ft, as in ©, by the essential 

bandwidth [—B/2,B/2] of the window g(t), and (3' is an enlargement of (3 by the support [— W/2,W/2] of the 

window g(t). Then, Q' = fl + B and ft = /? + W. 

Table U compares the number of samples necessary for a good approximation of time limited signals and of 

multipulse signals using these three methods. As can be seen from the table, the Gabor frame has two main 

advantages. The first is that it does not require strong decay of f(u>) for the reconstruction error to be bounded. 

Second, this approach can be used to efficiently sample multipulse signals with unknown pulse locations, as we 

will show in the next section. In this case we need approximately Afl'WNfi' 1 samples which is minimal with 

respect to the chosen approximation accuracy and frame redundancy. However, this amount increases slightly to 

the order of (9(4iV/i _1 \og(f3' / (4NW))£l'W) due to the utilization of CS algorithms in the recovery process. 
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Fourier 
series 


Shannon's 
interpolation 


Gabor series with 

S(g, a, 6)> ab = (j, 


number of 
samples 


» n"p 




ft fi'/J'/z- 1 


number of 
samples for 
multipulse 
signal with 
known pulse 
locations 


« n"/3 






number of 
samples for 
multipulse 
signal with 
unknown pulse 
locations 




ft Q"/3 


ft An'WNfi- 1 


approximation 
error 


® 


Q 


Theorem ImJl 



TABLE I: Comparison of three methods for approximating L2OR) functions that are time limited to [— /3/2,/3/2] 
and essentially bandlimited to [—17/2, £1/2]. The second and third lines refer to multipulse signals with N pulses, 
each of width no more than W. The methods are compared for the same accuracy of approximation. 

IV. Sampling of multipulse signals 

We now present a sampling scheme for functions from JAV(N, W, j3, f2) that reduces the number of channels 
in a Gabor sampling scheme and does not require knowledge of the pulse locations. 

A. Sampling System 

Our system, shown in Fig. Ufa), exploits the sparsity of multipulse signals in time. The signal f(t) enters JM 
channels simultaneously. In the (j, m)th channel, f(t) is multiplied by a mixing function q^ m [t), followed by an 
integrator. The design parameters are the number of channels JM and the mixing functions q^ m {t), < m < M— 1, 
< j < J — 1. The role of the mixing functions is to gather together all the information in f(t) over the entire 
interval [— /3/2, (3/2). Namely, f(t) is windowed with shifts of some compactly supported function, and all the 
windowed versions are summed with different weights. 

The functions qj, m (t) are constructed from the Gabor frame. Let Q(g, a, b) be a Gabor frame with window g(t) 
supported on the interval [—W/2,W/2], essentially bandlimited to [—B/2,B/2], and with sampling parameters 
a = /jW and b = 1/W for some < \i < 1. Then 

Qj,m(t) =Wj(t)s m (t) , (7) 
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Fig. 4: An efficient sampling system for multipulse signals (a), and an equivalent system using filters (b). The 
sampling step is r = WK. 



where 



^0 

/(*) = H d i' e " 

i=-L 



2-Kiblt 



k=-K 



with 



Sm(*) = ^ c mk9{t-ak), 

\n + B)w 



w 



2Wn 



1 and 



1 . 



(8) 



(9) 



Let K — 2Ko+l and L = 2Lq+1. The waveforms qj, m (t) are basically mixtures of K L channels e~ 2mblt g(t — ak), 
k = —Ko,.-.,Ko and I = — Lq, . . . ,Lq, of the Gabor sampling scheme, where the functions Wj(t) mix the 
frequency content of the signal, while s m (t) mix the temporal content of the signal. To specify qj im (t) completely, 
it remains to choose the coefficients dji and c m k defining the waveforms Wj(t) and s m (t), respectively. To do so, 
we first analyze the relation between the samples yj tTn and the signal f(t). 
Consider the (j, m)th channel: 



Vj,m — 



0/2 



0/2 



f(t)qj,m(t)dt 



A„ 



= X] d i l X c mk(f,M b iT ak g) 

I— — Lq h— — Kq 
L A 

1— — Lq k— — Ko 

The relation ( fTOb ties the known y m> i to the unknown Gabor coefficients z^.i of /(f) with respect to Q{g, a, b). This 
relation is key to the recovery of f(t). If we can recover Zk,i from the samples yj >m , then by Theorem IIII.ll we 



(10) 
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are able to recover f(t) almost perfectly. As can be seen from ( fTOb , the goal of the modulator qj, m (t) is to create 
mixtures of the unknown Gabor coefficients Zk,i- These mixtures, when chosen appropriately, will allow to recover 
Zk,i from a small number JM of samples by exploiting their sparsity and relying on ideas of CS. Note, that when 
using the basic Gabor scheme, each y j m is equal to one value of Zk,i, so that no combinations are obtained. When 
Zk.i are sparse, with unknown sparsity locations, we will need to acquire all their values using this approach. In 
contrast, obtaining mixtures of Zk.i, allows reduction in the number of samples. 

B. Signal Recovery 

For our purposes, it is convenient to write ( fTOb in matrix form as 

Y = DX T , with X = CZ. (11) 

Let the indices k = —Ko, ■ • • , Ko, I = —Lq, . . . , Lq, m = 0, . . . , M — 1 and j = 0, . . . , J — 1 be fixed throughout 
the exposition. Then, Y is a matrix of size J X M whose jmth element equals yj, m , and X is a matrix of size 
M x L with mZth element equal x m j. The unknown Gabor coefficients are gathered in the K x L matrix Z with 
columns Z[Z] = [z—k h ■ • • > z k ,i] T ■ The M x K matrix C contains the coefficients C m k+K — c mk, while the 
J x L matrix D contains the coefficients T)jj + L = dfl. The matrices C and D have to be chosen such that it is 
possible to retrieve Z from ( fTTT ). If J = L, M = K and D, C are identity matrices, then the system of Fig. Ua) 
reduces to standard Gabor sampling. 

From ( [Tol l it follows that the waveforms s m (t), respectively matrix C, mix the temporal content of the signal, 
while the waveforms Wj(t), respectively matrix D, mix the frequency content of the signal. The matrix C is used 
to reduce the number of channels. On the other hand, the purpose of D depends on which kind of signals are 
sampled. For general multipulse signals, the matrix D is only used to simplify hardware implementation, as we 
discuss below, but not to reduce sampling rate. Therefore, in general, we can choose D = I in this case. For 
multipulse signals that are additionally frequency sparse, we need D to allow recovery from lower rate samples, 
namely we can reduce the sampling rate by using appropriate mixtures with J < L, as shown in Section [V] 

We begin the discussion with general multipulse signals. When there is no frequency sparsity, we can choose 
J = L and D = I, reducing ( fTOb to X = CZ. In this case Wj(t) become pure modulations e ~ 27Ttb U- L o) _ Choosing 
J > L and D left invertible leads to a mixture of pure modulations, which can be easier to implement in hardware. 
This point is discussed in more detail in Section [VTTl Assuming D has full column rank, we can recover X from the 
samples Y by X = (DtY) T , where = (D ff D) _1 D H is the (Moore-Penrose) pseudoinverse of D. It remains 
to retrieve the unknown Gabor coefficients Zk,i from X = CZ. 

Recall from Fig. [3] that for every I, the column vectors Z[Z] of matrix Z have only \2fi~ 1 ~\N out of K nonzero 
entries, where the nonzero entries correspond to the pulse locations. In addition, all Z[Z] have nonzero entries in the 
same rows. The problem of recovering such a matrix Z is referred to in the CS literature as a multiple measurement 
vector (MMV) problem. Several algorithms have been developed that exploit this structure to recover Z efficiently 
from X in polynomial time when C has the RIP property of order 2[2//~ 1 ]iV, twice the number of nonzero rows 
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SUpp So(t) supp a^t + WK) supp S M -i(t + (M - l)WK) 

supp / 

m • • 

_^/2 £/ 2 1 i VTA- i i {M-1)WK ' * 

Fig. 5: Relation between the support of the filter s(t), which is the sum of the shifted supports of s m (t), and the 
support of the signal f(t). 



HZ], iTPJl . Ifl4l . ||26l . l27l . l28l . For example, a popular approach is by solving the convex problem 

min||Z|| 2 ,i subject to X = CZ , (12) 

where ||Z|| 2>1 = E£Lk (££- J*k,l| a ) 1/a - 

It is well known that Gaussian and Bernoulli random matrices, whose entries are drawn independently with equal 
probability, have the RIP of order S if M > cS\og(K/S), where c is a constant 11291 , ||30). For random partial 
Fourier matrices the respective condition is M > cS log 4 (K) QTI , ||32l . Therefore in our case, the number of 
samples in time has to be at least M > 2 [2/x _1 ] N \og(K/(2 \2/i~ 1 ~\ N j). 

C. Equivalent Representation 

For a fixed Gabor frame Q(g, a, b), the number of branches can be reduced to J if instead of JM modulations 
followed by an integrator, we perform J modulations followed by a filter s(t). Consider the system in Fig. HJb) 
with Wj(t) as in (0, r = WK, where K = 2Kq + 1, and Kq is as in (0, and the filter s(t) given by 

M-l 

s(t) = ^2 s m (t + WKm) . 
m=0 

Note, that for all to, s m (t) is compactly supported in time on [—W/2~fiWKo, W/2 + pWKo], and that its support 
contains the support [— /3/2,/3/2] of /(£). The shifted versions s m (t + WKm) have non-overlapping supports as 
the width of supp s m is smaller than the shift step WK 

W(l + 2fiK ) < W(l + 2K ) = WK . 

The support relation between the filter s(t) and the multipulse signal f(t) is depicted in Fig. [3] 
Under these assumptions, the output of the jth channel is 

yj, m = (wj(t)f(t) * s{-t))[WKm] 

M—l L a 

= djl{M-uf,T W K(m-n)S^)- 
n— l — — L(, 

The sum is nonzero only when m — n = 0, because otherwise the support of s n (t) shifted by WK{m — n) does 
not overlap the support of f(t), as depicted in Fig. Therefore it is sufficient to sample only at points t = WKm 
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for m = 0, . . . , M — 1, leading to 




^2 dji(M-bif,s m ) = E dji(f,M b is m ) 




Lo Ko 



E E 



l — — Lo k— — 



where Zk,i = (/, Mu T a k g) are the Gabor coefficients. Evidently, if the coefficients c/- TO used to build the blocks 
s m (t) of the filter s(i) are the same as coefficients used to create the waveforms qj. m (t), then the two systems are 
equivalent. 

D. Noisy Measurements 

Until now we considered signals that were exactly multipulse and noise free samples. A more realistic situation 
is when the measurements are noisy and/or the signal f(t) is not exactly multipulse, having some energy leaking 
outside the pulses. We now show that our sampling scheme is robust to bounded noise in both the signal and the 
samples. 

We say that a signal f(t) essentially bandlimited to [—£1/2,17/2], is essentially multipulse with N pulses each 
of width no more than W, if for some Sw < 1 there exists an f p £ AiViN, W,f3, O) such that 



We assume that the signals are time limited to the interval [—(3/2, (3/2], meaning that the energy leaks only between 
the pulses, and denote this class of signals by MVessiN, W, (3, il). 

Since the energy of / € AiVessiN, W, (3, f2) leaks beyond the support of the pulses, the column vectors Z[l], 
of the K x L matrix Z of dominant coefficients, defined in ( fTTI ). are no longer sparse. Nonetheless, Z can be well 
approximated by a sparse matrix Z s , which consists of S = \2/j,~ 1 ~\N rows of Z with largest £2 norm, and zeros 
otherwise, and is referred to as the best S— term approximation of Z. The existence of Z 5 is shown in Appendix 151 

Assuming now that the sampling system of Fig. Ufa) also has imperfections in the form of noise added to the 
samples, the input-output relation can be written as 



where X = CZ with Z a K x L matrix of Gabor coefficients and N is an J x M noise matrix. With D having 
full column rank, the relation dT3l > reduces to X = CZ + N, where N = D^N. A good S— term approximation of 
Z can be obtained by utilizing CS algorithms. Specifically, if C has RIP constant 62s < y2— 1 and N is bounded, 
then 



ll/-/p|| 2 < Ml/lb- 



Y = DX T + N 



(13) 



min||Z|| 2 ,i subject to ||CZ - X|| 2 < ||Nj| 2 



(14) 



has a unique S— sparse solution Z that obeys ll28l 



j|z-z|| 2 < d||z — z s ||2,i + c 2 ||isr||2 , 
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where C\ and C 2 are constants depending on <5 2 ,s. 

Finally, following a proof similar to that of Theorem lIII.il it can be shown that a function synthesized from Z is 
a good approximation of the original signal f(t): 

f- £ Zk,iM bl T akl \ < 

k— — K{) l — — Lo 



2 



< C (e B + e n )\\fh + Ci\\Z- Z*|| 2il + C 2 ||N|| 2 

where C = C^bllsollffllso. C 1 = C 0)6 ||7|| S( A, C 2 = C a , b || 7 || So C 2 and N = D+N. 

In particular, if Z is row sparse, as is the case for / e A4 V(N, W, (3, fi), then Z = Z s and the error of the 
approximation depends only on the noise added to the samples. When the signal is essentially multipulse, then the 
error bound depends on the decay of the coefficients. If that quantity is small, then a good approximation of f(t) 
is achieved by synthesizing a signal from the solution Z of (TT4l) . Note here, that the if the dual window 7 (t) is 
compactly supported, then a function reconstructed from the coefficients Z is multipulse. 

V. Time-Frequency sparse signals 

We now show that we can further reduce the sampling rate when sampling multipulse essentially multiband signals. 
We begin by giving a formal definition of such signals and describe the structure of their Gabor coefficients. 

A. Multipulse Essentially Multiband Signals 

We say that a signal / G A4V(N, W, j3, f2) is essentially multiband with S bands of width no more than ilw, 
if for some ew < 1 there exists a multiband function fi, with S bands, all of width no more then Qyy such that 

||/-/ 6 ||2<ew||./1|2. 

We denote the set of such signals by A4V(N, W, (3, S, ilw, An example are radar signals that are superpositions 
of a finite number of time-shifts and modulations of one pulse. If the generating pulse is well localized in frequency, 
then the signal is approximately sparse in the Gabor transform domain with respect to a window that decays fast 
in time and frequency. 

Let the Gabor frame be as in Section [IV] If the signal is known to be essentially multiband, then the nonzero 
row vectors Z[k] = [zk.~L , ■ ■ ■ , z k.L ] have only \(Cl\v + B)W]S out of L dominant entries, and the dominant 
entries correspond to the locations of the essential bands of f(t). Indeed, let [Bi,B 2 ] be any frequency band of 
/(*). Then B 2 - Bi < Q w and there are at most {(Qw + B)W] shifts of essential bandwidth [~B/2,B/2] of 
g(uj) by bl that overlap \B\, B 2 ]- This can be calculated from 

J2 ^ B_ _j_ Li 

2 2 W =^>L 1 -L 2 >(n w +B)W. 
Bi > # + # 

The coefficients for which the shift by bl of g(t) does not overlap any frequency band of f(t) are nonzero but small. 
Since there are altogether S bands present, each Z[fc] is \(flw + -B)W] S— dominant and all Z[fc] have dominant 
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entries on the same columns due to the structure of f(t). This implies that Z has at most \(Qw +B)W~\ S dominant 
columns, as shown in Fig. [3] 

The K x L matrix Zp with P = \(flw + B)W~\S nonzero columns corresponding to the P dominant columns 
of Z is referred to as the best P— column approximation of Z. Consequently, a result similar to that of Lemma IbTI 
holds with time and frequency interchanged. The sparsity in time combined with the sparsity of f(t) in frequency 
allows to further reduce the number of samples necessary for a good reconstruction. 

B. Signal Recovery 

To recover Z from the measurements in Fig. Ufa), the matrices C and D have to be chosen appropriately. If 
J > L and D is left invertible, then we are back to the situation of Section [TV] However, since Z is additionally 
almost sparse with respect to columns, we would like to reduce J. 

It is convenient to write the relation ( TTTb as 

Y T = CU , where U = ZD T . 

The matrix U inherits sparsity with respect to rows from the matrix Z, and therefore has only \2/i~ 1 ~\N out of 
K nonzero rows, which are precisely the nonzero rows of Z. When the matrix C has the RIP property of order 
\2(i~ 1 ~\N, then U can be efficiently recovered, for example by solving (TT2T > for a unique solution U subject to 
Y T = CU. 

Next, we use U to find a unique \(fl\y + B)W~\ S— sparse approximation of Z. Let S be the set of indices of 
nonzero rows of U and U 5 the matrix built from those rows of U indexed by S. If the matrix D has RIP constant 
S2P < V2 — 1 with P = \(Qw + -B) W] S, then there exists a unique P~ sparse solution V of 

iriin|| V|| 2 ,i subject to (U 5 ) T = DV . 

Let Z be a K x L matrix whose \2fj,~ 1 ~\N rows indexed by S equal to V T , and the remaining entries equal to 
zero. It then follows that Z is proportional to the best P— column approximation of Z in the following sense |28| 

||Z-Z|| 2 < Ci\\(Z-Z v ) T \\ 2 ,i, 

where C\ is a constant depending on <52p. The requirement on D translates to J > 2 \(£lw + B)W~\ S. As opposed 
to purely multipulse signals, where it suffices to take D = I, this choice is not possible here, since it does not 
satisfy RIP. 

The resulting matrix Z is a [2/i -1 ] N row sparse and [(0,w + B)W~\S column sparse approximation of Z. It 
is important to note, that the solution to the first MMV problem, Y T = CU, recovers U exactly, since U is row 
sparse, while the solution to the second MMV problem, U T = DZ T , returns a column sparse matrix Z that is an 
approximation of Z, which itself is not strictly column sparse. 
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Finally, the function reconstructed from the coefficients Zk,i is a good approximation of the input signal f(t): 

K L 
f - } J ^ Zk,lM b iT ak J ^ < 
k— — Kq I — — Zjq 

< C (e B + e n )\\fh + Ci\\(Z - Z V ) T \\ 2A 

where Cq — C\ (JMIsolMlso an d @i = Cq;,6||tIIsoCi- The proof is analogous to the proof of the error estimate in 
Theorem IIH. II with appropriate adjustments. 

In the case of known positions of the pulses and bands, the minimal sampling rate for the desired accuracy of 
the approximation and a given frame is when M — \2fi~ 1 ~\ N and J = \(ftw + B)W~\ S. In the blind setting, when 
the locations of the pulses and the bands are unknown, the sampling rate increases by a factor of four (a factor of 
two in each domain), with M > 2\2/j,~ 1 ~\N and J > 2[(f2yK + B)W~\S required for obtaining a unique solution. 
Therefore, for signals from the set .MT^iV, W, f3, S, fl\y, the number of samples M with respect to time is the 
same as for signals from AiViN, W, (3, il), while J, the number of samples with respect to frequency, is reduced 
from J > L to J < L. The overall number of samples is MJ « 8n,' w W N S fi^ 1 , where Vt' w = Vlw + B. 

VI. Related work 

Recently, the ideas of CS have been extended to allow for sub-Nyquist sampling of analog signals [6 1, [8 1, [9 1, [17], 
ifTSl , [28 1, [33 1, l34l . These works follow the Xampling paradigm, which provides a framework for incorporating 
and exploiting structure in analog signals without the need for discritization |fT31 , fl6l . Two of these sub-Nyquist 
solutions are closely related to our scheme: the first is a sub-Nyquist sampling architecture for multiband signals 
introduced in IfTTl , while the second is a sampling system for multipulse signals with known pulse shape introduced 
in [8|. We show, that by choosing different waveforms g j m (t), the systems of IfTTl and ||8) are special cases of the 
system in Fig. [4] 

A. The Modulated Wideband Converter 

The concept of using modulation waveforms is based on ideas presented in IfTTl for a multiband model, which 
is Fourier dual to ours: the signals in [17| are assumed to be sparse in frequency, while multipulse signals are 
sparse in time. More specifically, [17| considers multiband signals whose Fourier transform is concentrated on TV 
frequency bands, and the width of each band is no greater than B. The locations of the bands are unknown in 
advance. A low rate sampling scheme, called the modulated wideband converter (MWC), allowing recovery of such 
signals at the rate of ANB was proposed in IfTTl ; a hardware prototype appears in |fT6l . This scheme consists of 
parallel channels where in each channel the input is modulated with a periodic waveform followed by a low-pass 
filter and low-rate uniform sampling. The main idea is that in each channel the spectrum of the signal is scrambled, 
such that a portion of the energy of all bands appears at baseband. Therefore, the input to the sampler contains a 
mixture of all the bands. Mixing of the frequency bands in IfTTl is analogous to mixing the Gabor coefficients in 
our scheme. 



December 21, 2011 



DRAFT 



20 



The MWC is equivalent to the system of Fig. |4jb) where the waveforms Wj(t) are B— periodic and the filter 
s(—t) is an ideal rectangular low pass filter, whose bandwidth is [— B/2, 5/2]. The samples are taken at points 
t = m/B, m E Z. The output of the MWC system is then a weighted sum of Gabor coefficients with respect to 
a frame Q(g r ,a,b) where g r (t) is a sine function that is bandlimited to [— B/2, B/2] and a = l/B, b = B. Thus 
the samples can be written as 

Lo 
l=-L 

where Lq = + B)/ (2B)] — 1. With this frame, for each m € Z the number of nonzero Gabor coefficients z m j 
equals at most 2B, as at most two shifts of [— B/2, B/2] by bl = Bl overlap one band of the signal. Therefore, 
the number of channels J is proportional to the number N of frequency bands in the signal, and equals J > ANB. 

The MWC is an ideal system, in the sense that it uses ideal low pass filters, which in practice are difficult to 
build, and that the reconstruction process uses infinitely many samples. Using the Gabor approach we can generalize 
the MWC to other, not necessarily ideal filters. Furthermore, the reconstruction error can be computed when only 
a finite number of samples is available by using Theorem Mill with time and frequency interchanged. 

The MWC can be easily extended to other, more redundant frames, with a cost of increased number of channels 
J. Let Q(g,a,b) be a collection of Gabor frames with windows g(t) bandlimited to [— B/2, B/2] and constants 
a = l/B and b = pB, for some p e (0,1). For any frame from Q{g,a,b), the MWC parameters change to 
Lq = + B) /(2Bp)~\ — 1, the waveforms Wj(t) have to be Bp— periodic and the filter s(t) — g(t). Also, the 
sparsity of Gabor coefficients in frequency reduces to \2p~ 1 ~\N, as |~2/i _1 ] shifts of [— B/2, B/2] by bl = Bpl 
overlap one band of f(t). The MWC system associated to this frame has to have J > 2|~2 / u _1 ] NBp -1 channels. 
This is an increase in the number of channels by a factor of p, . However, this increase can reduce the number 
of time samples necessary for achieving the same reconstruction error as with Q(g r , 1/B,B). 

B. Multipulse Signals with Known Pulse Shape 

Another related signal model is that of multipulse signals with known pulse shapes 0, El, J9): 

s 

f(t)=J2^h(t-t s ) (15) 

where h(t) is known and f(t) is supported on [—[3/2, (3/2]. This problem reduces to finding the amplitudes a s and 
time delays t s . As shown in [6| the time-delays can be estimated using nonlinear techniques e.g. the annihilating 
filter method as long as the number of measurements L satisfies L > 2S and the time-delays are distinct. Once 
the time delays are known, the amplitudes can be found via a least squares approach. The number of channels is 
motivated by the number of unknown parameters (<r s ,t s ) which equals 2S. 

The Fourier coefficients can be determined from the samples of f(t) using a scheme similar to that of Fig. 2Ja) 
with J > L channels. The J modulating waveforms being q^ m (t) = Wj(t) with b = 1/(3, and all s m (t) set to one. 
In this case, the input-output relation becomes 

y = Df, 
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where y is a vector of length J, D is a matrix of size J x L and f is a vector of Fourier coefficients f(l//3) of 
f(t) of length L. If Wj(t) are designed so that D is left invertible, then f = D^y. We note here, that the system of 
1 8 1 is inefficient for our signal model, since it reduces to the Fourier series method, which does not take sparsity 
in time into account. However, by choosing an appropriate Gabor frame and waveforms qj >m (t), the same scheme 
of Fig. a) can be used both to sample signals from the set M.V(N, W, (3, il) as well as that of the form (TT3T >. as 
shown in the following proposition. 

Proposition VI.l. Let Q{g,a,b) be a Gabor frame such that X)fcezff(^ — ~ 1 almost everywhere, and the 
waveforms qj_ m (t) in the sampling scheme of Fig.^aj are such that the matrix D, is left invertible and the matrix 
C has RIP constant 5 < y2 — 1 of order 2\2fi~ 1 ~\ N with one row of ones. Then this sampling scheme can be used 
to sample multipulse signals of the form ( 1751 ) supported on [— (3/2, /3/2]. The time-delays and amplitudes of f(t) 
can be retrieved from J samples as long as t s £ [— W/2, W/2], L > 2S and the known pulse h(t) in ( 1751 ) satisfies 
h(bl) ^Ofor \£\ < L Q . 

The proof is straightforward. 

Example of Gabor windows g(t) that are well localized in time and frequency and form a partition of unity, e.g. 
J2k£Z,d(t ~ = 1' aie me ra i se d cosine window, or B-splines of positive orders ||23l . An example of a matrix 
C with a row of ones is a partial Fourier matrix which is known for its good CS properties |2TI . 

To conclude, we have seen that the same hardware can be used to sample signals with known pulses and those 
from AiViN, W, (3, il). The difference is in the number of branches used and the processing stage. 

VII. Waveform Design 

Hardware implementation of our scheme reduces to implementing the waveforms qj, m {t). The mixing functions 
qj m (t) are a product of Wj(t) and s m (t) defined in ((HJ. The functions s m (t) are pulse sequence modulations, 
where the sequences are generated to form a valid CS matrix. An example is a matrix whose entries are ±1 drawn 
independently and with equal probability. 

One method to create the waveforms Wj(t) is to low-pass filter 1/6— periodic waveforms. More precisely, let 

l£Z i=0 V ' 

where p(t) is some pulse shape such that w(bl) ^ for \£\ < Lq, and aj[i] is a length— I sequence. Since Wj(t) 
is 1/6 periodic it can be expressed as 

for some coefficients dji. We then filter Wj(t) by a filter u(t) with frequency response u(lu), designed so that 

{1 UJ = bl, |Z| < L 

u = bl, \l\ > L Q , 

arbitrary elsewhere 
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to form the waveforms Wj(t) — Y^t=-L djie~ 27rMt with coefficients dji — dji -u(bl). The shaping filter frequency 
response, u(uj), is designed to transfer only the coefficients with index £ — —Lq, . . . ,Lq, suppressing all other 
coefficients. 

For the matrix D, built from the coefficients dji, to be left invertible a necessary condition is that J > I > L 
and the sequences ctj[i] are chosen such that the matrix A, whose jith element is OLj[i], has full column rank [8|. 
For example, if J = I = L, then the rows of A can be created from cyclic shifts of one basic sequence. On the 
other hand, for a matrix D to be a valid CS matrix, meaning to have RIP property with high probability, the values 



a,- i 



±1 are chosen independently with equal probability and I > L > J |fT71 . 



w(t) 



One example of a pulse modulation scheme is when J = I = L, and 

1 *e [o,&] 

The frequency response of this pulse is given by 

w{u) = ±e-ff ■ sine (£) , 

so that w(bl) 7^ for \£\ < L . In addition we choose cnj\i\ as sequences of ±ls, created from cyclic shifts of one 
basic sequence, in a way that yields an invertible matrix A. Such rectangular pulses with alternating signs can be 
easily implemented in hardware [16|. 



VIII. Gabor windows 

The sampling scheme presented in this paper is based on Gabor frames. We recall here some methods to construct 
Gabor frames with well localized windows for a chosen redundancy /i based on results from [22] and [23 1. 

Let fj, > 4. A window g(t) that is supported on [— W/2, W/2] and forms a frame with a = [iW and b = 1/W 
can be constructed from an everywhere increasing function h(t) such that h(t) = for t < 0, and h(t) = 1 for 
t > 1 by 

' 0, 



9(t) = { 



1/2 



l _ h( tJW-X I 2 



1/2 







t€[-f 


wx 

2 


\t\<^ 






W] 
2 J ' 


*>¥. 





0, 

where A = 2/i— 1, 1221 . The function g(t) is non-negative, has the desired support and equals 1 on [— W\/2 1 WX/2]. 
If h(t) is taken to be 2k continuously differentiable, than g(t) is k times continuously differentiable, which implies 
that g(i>j) decays like o(|w| ). The points t = ±WA/2, where g(t) becomes constant, have been chosen so that 
their distance to the furthest edge of suppg is exactly /j,W. The frame bounds of such a constructed frame equal 
Ai = A 2 = 1 E2, since 

5>(i+ k»W)\ 2 = l. 

fcGZ 
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As an example, let h(t) = sin(7rf/2) on [0, 1] and /x = i. Then 

M / \t\> W/2, 

[ cos{nt/W) \t\ < W/2. 

An alternative construction for /j, > 1/2 was developed in [35], [36|. The method results in spline type windows 
g(t) of any order of smoothness that satisfy the partition of unity criterion. The constructions are made by counting 
the number of constraints (in the Ron-Shen duality condition |[3"71 , and on the points where continuity /differentiability 
is required) and then searching for polynomials on [— 1, 0] and on [0, 1] of a matching degree. One example is g(t) 
supported on [—2/3, 2/3] and given by 

{2 + 3* t G [—2/3, —1/3] , 
1 |*[ < 1/3, 

2-3* te [1/3,2/3], 

that forms a frame with a = 1 and b = 3/4. It forms a partition of unity with a shift parameter a = 1, J^kezdtt ~ 
k) = 1. The dual window is also supported on [—2/3, 2/3] and is given by 

{-18* 2 -15t-2 [-2/3,-1/3], 
1 |*| < 1/3, 

-18* 2 + 15t-2 te [1/3,2/3]. 

Applying dilation by (/iW) -1 , with /i = 3/4, both to g(t) and 7(*) we obtain a dual pair of windows g(t) and 
lit) 

g(t) = g(t/(pW)) 7 (i) = 7W(MW0) 

that are supported on [— W/2, W/2], and such that Q(g.jiW,l/W) forms a frame with frame bounds A\ = 1/2 
and A2 = 1. Moreover, g(t) forms a partition of unity with shift parameter a = [iW , J2kez 9(t ~ iWVk) = 1. 
Well known, compactly supported Gabor windows are the B-splines. Let -Bjv(*) be a spline of order N, 

Bi(t) = Xi/ 2 (*) , B N+1 (t) = (B N * Bi)(t) . 

Then i?Ar(t) is supported on [—AT/2, AT/2] and forms a partition of unity with shift parameter a = 1. To generate 
a Gabor frame from i3jv(*) with a window supported on [—W/2, W/2] and lattice parameters a = fiW, b = l/W, 
such that the window forms a partition of unity with shift /j,W, we need to choose /i — l/N l38l . Then g(t) = 
B]y(tN/W) is supported on the desired interval and decays like (1 + |cj|) _Ar_c in the frequency domain. Note that 
fi decreases as the order N of smoothness of the B-spline is increased. Thus smoother windows can be obtained 
only at the cost of a smaller /i. However, already for N = 3 we get good concentration properties of g(t). The 
dual can be computed by inverting the Gabor frame operator, or by using the method of l23l . 

IX. Simulations 

We now present some numerical experiments illustrating the recovery of multipulse signals. 
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Fig. 6: (a) Decay of the reconstruction error with increased number of samples for multipulse signals with N = 1,3,5 
pulses, pulse width W = 0.18ms and frame redundancy fi = 0.5. (b) Comparison of performance for signals with 
N = 3, W — 0.24ms using different frames with /i = 0.3, \i = 0.5 and /i = 0.75. (c) Comparison of the 
relative error with respect to the number of channels J for multipulse essentially bandlimited signals with N = 3, 
W = 0.18ms, S = 2 with respect to different frames. 



We tested our sampling scheme using Monte Carlo simulations averaged over 500 trials on a range of multipulse 
signals of duration j3 = 22ms. The pulses making up the signals were chosen at random from a set of five different 
pulses: cosine, Gaussian, B-spline of order 3 and 5, and rectangular pulse. The locations of the pulses were also 
chosen at random. We varied the number of pulses N = 1,3,5, the maximal width W = 0.18,0.24ms of the 
pulse, and the redundancy /i = 0.3, 0.5, 0.75 of the frame. Throughout the experiments we chose D = I and C as 
a Bernoulli random matrix. We measured the relative error ||/ — /||/||/||2- For redundancy fi — 0.3 we chose a 
Gabor frame with window being B-spline of order three, for /i = 0.5 the window was a cosine, and for /i = 0.75 
we chose the truncated Gaussian. 

Fig-I2f a ) depicts the decrease in the reconstruction error with increased number of samples M for different values 
of N and W = 0.18ms. We used a tight Gabor frame with a cosine window and redundancy fj, = 0.5. The M is 
dictated by the the number of pulses and frame redundancy, and it has to be at least 2[2/i~ 1 ]iV. Meaning, that for 
multipulse signals with N = 1 pulses, M > 8, for N = 3 we have M > 24, and for N = 5 it has to be M > 40. 
As expected, the sparser the signal, the less samples are needed for a good reconstruction. The number of samples 
in time can be significantly reduced if sparsity is taken into account. Without any knowledge on the sparsity we 
would have to take K = 241 time samples for signals with N = 5 pulses, and with Lq = 5 that would result in 
the reconstruction error of 0.05. However, when sparsity is taken into account, already M — 40 samples suffice to 
achieve the same reconstruction error. Therefore reducing the number of samples by a factor of six. When N = 3, 
to achieve reconstruction error of 0.05 we need M = 28 samples in time, and M = 12 for signals with N = 1. 

In Fig. |6jb) we considered the influence of the Gabor frame on the reconstruction error and the number of 
samples M involved. We tested the system for signals with N — 3 pulses of width no more than W — 0.24ms and 
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Lq = 5. The least number of samples K = 121 is achieved with /i = 0.75 and at the same time with M = 25 we 
achieve a good reconstruction. The value of M necessary for a good reconstruction increases with the increase of 
redundancy. Without knowing the sparsity structure of the signal in time, we would have to take K = 271 samples 
for /i = 0.3, and K = 120 for fi — 0.5. When sparsity is exploited, we can reduce that number to M — 45 and 
M = 27, respectively. 

We then examined the performance of our sampling scheme on signals comprising three pulses of width no 
more than W = 0.18ms, that are additionally essentially multiband with two bands. Fig. |6fc) depicts the decay of 
reconstruction error with the increase of J for two different frames: one is a tight frame with cosine window and 
redundancy /i = 0.5 and second, a frame with Gaussian window and redundancy fi = 0.75. The sampling system 
was tested with the matrix D being the random Fourier matrix and C a Bernoulli random matrix. For example, 
when a frame is of redundancy /i = 0.5 and no sparsity is taken into account then we need K = 241 and L = 91 
samples to achieve a reconstruction error of 0.07. On the other hand, with sparsity being exploited we can use 
only M — 35 and J — 40 for a similar reconstruction quality, resulting in a twelvefold reduction in the number of 
samples. 

X. Conclusions 

We presented an efficient sampling scheme for multipulse signals, which is designed independently of the time 
support of the input signal. Our system allows to sample multipulse signals at the minimal rate, far below Nyquist, 
without any knowledge of the pulse shapes or its locations. The scheme fits into the broad context of Xampling - 
a recent sub-Nyquist sampling paradigm for analog signals. Our architecture relies on Gabor frames which lead to 
sparse expansions of multipulse signals, and consists of modulating the signal with several waveforms followed by 
integration. We showed that the Gabor coefficients, necessary for reconstruction, can be recovered from the samples 
of the system by utilizing CS techniques. The number of necessary samples depends on the desired accuracy of the 
approximation, essential bandwidth of the signal, and redundancy factor fi related to the Gabor frame, and equals 
A^I'NW/j,^ 1 . The sampling rate can be further reduced if the signal is additionally sparse in frequency. We also 
showed that the proposed sampling and recovery technique is stable with respect to noise and mismodeling. 

Appendix A 
Proof of Theorem IIII.1I 

The proof is rooted in that of Theorem 3.6.15 in [24 1 with appropriate adjustments. Since Q(g,a,b) is a Gabor 
frame, f(t) admits a decomposition 

Kg 

/ = zk,iM b i T ak 7 . 

k=-K zez 

Let cb > 0. The bandlimited So functions are dense in So, therefore, there exists g c € So bandlimited to some 
[-5/2,5/2], such that 

ll.g-.gdl s < es||.g||s - 
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Since f(t) is an essentially bandlimited function, there exists a function f c (t) bandlimited to [— f2/2,Q/2], such 
that 

||/-/c||2<eo||/|| 2 . 

Consequently, \zk,i\ — \ (f c , M- a k Tug c )\ ^ only for those I such that supp/ c n (supp ^ + bl) ^ 0, that is 

[-0/2, il/2] n [bl - B/2, bl + 5/2] ^ . 

The fact that f c (t) and g c (t) are bandlimited implies that there are only a finite number of values I for which 
Vg c f c (ak, bl) ^0. Let Lo be the smallest integer such that \V gc f c {ak 1 bl) \ — for |2| > Lo. The exact value of Lo 
can be calculated as 



Lo 



n + B 

2b 



1 . 



Define a sequence dk,i as 



<4j = 



Zk.i , \k\ < K , \l\ > L Q 



, else. 

Then |d fc ,i| < \V g - g J(ak, bl) + V 9c (f - f c )(ak, bl)\ for all k, I G Z, and 

/- £ S Z M M «T afe7 
k— — KQ 1=—Lq 

fcez zez 

<Ca,6||7l|So(IIV flc /|k + ||V ffe (/-/ c )|| fe ) 

<c a 2 b || 7 |U ||g||5o(^ + ^)ll/ll2 

where we first used the boundedness of the analysis operator related to g(t) and then the synthesis operator related 
to j(t) whenever g and 7 are in So- 



Appendix B 

S*— TERM APPROXIMATION OF Z 

We show here the existence of an S— term approximation of Z. 

Lemma B.l. Let f G A4V ess (N, W, (3, f2) be 8w —essentially multipulse and Q(g,a,b) be a Gabor frame with g 
compactly supported on [-W/2, W/2] and a = pW, b = 1/W for some < // < 1. Then there exists a subset S 
of {—Ko, . . . , Kq} such that 



||Z-Z*|| 2 ,i <6wVKC a , b \\g\\ So \\fh, 
where Z s consists of rows of Z indexed by S, K = 2K + 1 and ||Z|| 2 ,i = J2k=-Ko^^i J =-L \ Zk ^ 2 ^ 1 ^ 2 - 

Proof: Let f p G M.V(N, W, (3, ft) be a multipulse 5w— approximation of /. Then V g f p (ak,bl) = for all 
\k\ > Kq, and the column vectors [V g f p (— aK , bl), . . . , V g f p (aK , bl)] T , \l\ < L , are all jointly sparse with 
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\2[i~ 1 ]N nonzero coefficients. Let S denote the index set of nonzero coefficients. For \£\ < L , let Z s [l] be 
vectors with coefficients z^ l defined by 

Zk,i k e S 
k £ S. 

Then Z s [l] is the best |"2/x _1 ]./V-term approximation of Z[Z], for each 1. Note that \z k .i-zfj\ < \V g (f- f p )(ak, bl)\ 
for all k and £, so that 

iiz-z 5 n 2 , 1= ( jr \z kA -ziA 

k=-K a \l = -L / 
K / L \ 1/2 

^ E E i^(/-/ P )(«*.w)i 2 

fe=-if \l = -L / 
K 

< E \\V g (f-fp)(ak,-)\U,<VK\\V g (f-f p )\\ t3 

k=-K 

<VKC a M\g\\s \\f-f P h<5 w VKC a M\g\\s a \\fh, 

completing the proof. ■ 
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